/*-
 * Copyright (c) 2013-2014 Shang-Rong Cai. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHORS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 *
 * $ pregrid.c    Rev. 0.2.2    2014-05-30T11:03+08:00 $
 */

#include <stdlib.h>
#include <stdio.h>

#include <pregrid.h>

int cidk_init(size_t rz, size_t cz, size_t kz, uint8_t flg, cidk_t *bk)
{
    int ierr = ERRNO_UNDEFINED;
    size_t ptrz;

    ptrz = 0;

    bk->r_sz = rz;
    bk->c_sz = cz;
    bk->e_sz = 0;
    bk->k_sz = kz;

    bk->flag = flg;

    switch (bk->flag)
    {
        case ROW_WISE:
            ptrz = bk->r_sz;
            break;
        case COL_WISE:
            ptrz = bk->c_sz;
            break;
        case RDR_NONE:
            ierr = ERRNO_SUPPORT_NONE;
            goto EXIT_CIDK_INIT;
            break;
        default:
            ierr = ERRNO_UNKNOWN_TYPE;
            goto EXIT_CIDK_INIT;
            break;

    }

    bk->sptr = (sidxpair_t*) calloc(ptrz, sizeof(sidxpair_t));
    bk->item = (cidbl_t*) calloc(bk->k_sz, sizeof(cidbl_t));

    EXIT_CIDK_INIT:

    return ierr;
}


int cidk_dest(cidk_t *bk)
{
    int ierr = ERRNO_UNDEFINED;

    bk->r_sz = 0;
    bk->c_sz = 0;
    bk->e_sz = 0;
    bk->k_sz = 0;

    bk->flag = RDR_NONE;

    free(bk->item);

    return ierr;
}


int cidk_stdout(cidk_t *bk)
{
    int ierr = ERRNO_UNDEFINED;

    size_t zi;
    sidx_t si;

    size_t cmpidxz;  /* Compressed index size */
    size_t cnt;

    switch (bk->flag)
    {
        case ROW_WISE:
            cmpidxz = bk->r_sz;
            break;
        case COL_WISE:
            cmpidxz = bk->c_sz;
            break;
        case RDR_NONE:
            ierr = ERRNO_SUPPORT_NONE;
            goto EXIT_CIDK_STDOUT;
        default:
            ierr = ERRNO_UNKNOWN_TYPE;
            goto EXIT_CIDK_STDOUT;
            break;
    }

    cnt = 0;

    if (bk->flag == ROW_WISE)
    {
        for (zi = 0; zi < cmpidxz; zi++)
        {
            for (si = bk->sptr[zi].head; si < bk->sptr[zi].tail; si++)
            {
                printf("%6zd: (%6zd, %6d) % .14le\n", cnt, zi, bk->item[si].index, bk->item[si].value);
                cnt++;
            }
        }
    }
    else
    {
        for (zi = 0; zi < cmpidxz; zi++)
        {
            for (si = bk->sptr[zi].head; si < bk->sptr[zi].tail; si++)
            {
                printf("%6zd: (%6d, %6zd) % .14le\n", cnt, bk->item[si].index, zi, bk->item[si].value);
                cnt++;
            }
        }
    }

    EXIT_CIDK_STDOUT:

    return ierr;
}


/* Matrix Vector Multiplication y=Ax */
int cidk_vector_multi(const cidk_t *cidk, const double *vec_x, double *vec_y)
{
    int ierr = ERRNO_UNDEFINED;

    size_t zi;
    sidx_t si;

    for (zi = 0; zi < cidk->r_sz; zi++)
    {
        vec_y[zi] = 0;

        for (si = cidk->sptr[zi].head; si < cidk->sptr[zi].tail; si++)
        {
            vec_y[zi] += cidk->item[si].value * vec_x[cidk->item[si].index];
        }
    }

    return ierr;
}


/* Two Dimension Domain Cells */
int shpfn_std_tria(int id, char pr, double vr, double vs, double *vl)
{
    int ierr = ERRNO_UNDEFINED;

    switch (id)
    {
        case 0:
            switch (pr)
            {
                case 'o':
                    *vl = 1.0 - vr - vs;
                    break;
                case 'r':
                    *vl = -1.0;
                    break;
                case 's':
                    *vl = -1.0;
                    break;
                default:
                    ierr = ERRNO_UNKNOWN_STYLE_PARTIAL;
                    break;
            }
            break;
        case 1:
            switch (pr)
            {
                case 'o':
                    *vl = vr;
                    break;
                case 'r':
                    *vl = 1.0;
                    break;
                case 's':
                    *vl = 0;
                    break;
                default:
                    ierr = ERRNO_UNKNOWN_STYLE_PARTIAL;
                    break;
            }
            break;
        case 2:
            switch (pr)
            {
                case 'o':
                    *vl = vs;
                    break;
                case 'r':
                    *vl = 0;
                    break;
                case 's':
                    *vl = 1.0;
                    break;
                default:
                    ierr = ERRNO_UNKNOWN_STYLE_PARTIAL;
                    break;
            }
            break;
        default:
            ierr = ERRNO_OUT_OF_RANGE_SHAPE_ID;
            break;
    }

    return ierr;
}


int refelm_tria_q1(refelm_t *refelm, void *opfn)
{
    int ierr = ERRNO_UNDEFINED;
    void (*opcode)(int, char, double, double, double*);

    int ii, ij;  /* int based loop counter or index */

    double qp_x;  /* Integration points x component */
    double qp_y;  /* Integration points y component */

    int cur_index;  /* Shift index */

    /* Setup Coordinates And Weights of Quadrature Points */
    qp_x = 1.0 / 3.0;  /* Set integration point 0 for x component */
    qp_y = 1.0 / 3.0;  /* Set integration point 0 for y component */
    refelm->qdtwgt[0] = 1.0 / 2.0;  /* Set quadrature weights */

    /* Set Reference Shape Functions */
    opcode = opfn;

    /* Setup Reference Values */
    for (ii = 0; ii < refelm->bssfnz; ii++)
    {
        cur_index = ii * refelm->qdtptz;

    #ifdef DEBUG_INFO
        printf("Local Node %d:\n", ii);
    #endif

        for (ij = 0; ij < refelm->qdtptz; ij++)
        {
            opcode(ii, 'o', qp_x, qp_y, &refelm->refval[cur_index + ij]);
            opcode(ii, 'r', qp_x, qp_y, &refelm->refpvr[cur_index + ij]);
            opcode(ii, 's', qp_x, qp_y, &refelm->refpvs[cur_index + ij]);

        #ifdef DEBUG_INFO
            printf("             val: %lf\n", refelm->refval[cur_index + ij]);
            printf("             pvr: %lf\n", refelm->refpvr[cur_index + ij]);
            printf("             pvs: %lf\n", refelm->refpvs[cur_index + ij]);
        #endif
        }
    }

    /* Finalize */
    EXIT_REFELM_TRIA:

    return ierr;
}


int refelm_tria_q3(refelm_t *refelm, void *opfn)
{
    int ierr = ERRNO_UNDEFINED;
    void (*opcode)(int, char, double, double, double*);

    int ii, ij;  /* int based loop counter or index */

    double *qp_x;  /* Integration points x component */
    double *qp_y;  /* Integration points y component */

    int cur_index;  /* Shift index */

    /* Initialize */
    qp_x = (double*) calloc((size_t) refelm->qdtptz, sizeof(double));
    qp_y = (double*) calloc((size_t) refelm->qdtptz, sizeof(double));

    /* Setup Coordinates And Weights of Quadrature Points */
    qp_x[0] = 0;  /* Set integration point 0 for x component */
    qp_y[0] = 0;  /* Set integration point 0 for y component */
    refelm->qdtwgt[0] = 1.0 / 6.0;  /* Set quadrature weights */

    qp_x[1] = 1.0;  /* Set integration point 1 for x component */
    qp_y[1] = 0;  /* Set integration point 1 for y component */
    refelm->qdtwgt[0] = 1.0 / 6.0;  /* Set quadrature weights */

    qp_x[2] = 0;  /* Set integration point 2 for x component */
    qp_y[2] = 1.0;  /* Set integration point 2 for y component */
    refelm->qdtwgt[0] = 1.0 / 6.0;  /* Set quadrature weights */

    /* Set Reference Shape Functions */
    opcode = opfn;

    /* Setup Reference Values */
    for (ii = 0; ii < refelm->bssfnz; ii++)
    {
        cur_index = ii * refelm->qdtptz;

    //#ifdef DEBUG_INFO
        printf("Local Node %d:\n", ii);
    //#endif

        for (ij = 0; ij < refelm->qdtptz; ij++)
        {
            opcode(ii, 'o', qp_x[ij], qp_y[ij], &refelm->refval[cur_index + ij]);
            opcode(ii, 'r', qp_x[ij], qp_y[ij], &refelm->refpvr[cur_index + ij]);
            opcode(ii, 's', qp_x[ij], qp_y[ij], &refelm->refpvs[cur_index + ij]);

        //#ifdef DEBUG_INFO
            printf("             val: %lf\n", refelm->refval[cur_index + ij]);
            printf("             pvr: %lf\n", refelm->refpvr[cur_index + ij]);
            printf("             pvs: %lf\n", refelm->refpvs[cur_index + ij]);
        //#endif
        }
    }

    /* Finalize */
    EXIT_REFELM_TRIA:
    free(qp_x);
    free(qp_y);

    return ierr;
}


int refelm_triangle_init(refelm_t *refelm)
{
    int ierr = ERRNO_UNDEFINED;
    size_t ref_vz;

    refelm->dimz = 2;  /* 2D */

    refelm->qdtptz = 1;  /* Quadrature point */
    refelm->qdtwgt = (double*) calloc((size_t) refelm->qdtptz, sizeof(double));

    refelm->bssfnz = 3;  /* Triangle */

    ref_vz = (size_t) (refelm->bssfnz * refelm->qdtptz);

    refelm->refval = (double*) calloc(ref_vz, sizeof(double));
    refelm->refpvr = (double*) calloc(ref_vz, sizeof(double));
    refelm->refpvs = (double*) calloc(ref_vz, sizeof(double));

    ierr = refelm_tria_q1(refelm, (void*) shpfn_std_tria);

    return ierr;
}


int refelm_triangle_dest(refelm_t *refelm)
{
    int ierr = ERRNO_UNDEFINED;

    free(refelm->qdtwgt);

    free(refelm->refpvr);
    free(refelm->refpvs);

    return ierr;
}
